#include "mvn.h"

double InvErrFunc(double x)
{
	double result;
	try
	{
		result = boost::math::erf_inv(x);
	}
	//Check for input values: 1.0 and -1.0
	//In these cases, an exception will be thrown
	catch(exception &)
	{
		if(x>0)
		{
			result = 1e+300;
		}
		else
		{
			result = -1e+300;
		}
	}
	return result;
}


Matrix<double> CholeskyFactor(const Matrix<double> &mat)
{
	size_t n=mat.getSize();
	Matrix<double> c(n);
	for(size_t i=0 ; i<n; i++){
		for(size_t k=0; k<(i+1); k++){
			double sum=0.0;
			for(size_t j=0; j<k; j++){
				sum+=c[i][j]*c[k][j];
			}
			c[i][k]=(i==k)?sqrt(mat[i][i]-sum):(1.0/c[k][k]*(mat[i][k]-sum));
		}
	}
	return c;
}